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Abstract — Recent  work  on  block-based  compression  for  low 
bit-rate  coding  has  shown  that  employing  a  block  coder  within  a 
sampling  scheme  where  the  image  is  downsampled  prior  to  coding 
(and  upsampled  after  the  decoding  stage)  results  in  superior  per¬ 
formance  compared  to  standard  block  coding.  In  this  paper,  we 
explore  the  use  of  optimal  decimation  and  interpolation  filters  in 
this  coding  scheme.  We  show  that  the  problem  of  finding  optimal 
filters  for  a  general,  unknown,  “black-box”  coder  can  be  written 
as  a  separable  least  squares  problem  in  two  sets  of  variables.  We 
then  elegantly  solve  this  optimization  problem  using  the  Variable 
Projection  method.  The  experimental  results  presented  clearly 
exhibit  a  significant  improvement  over  existing  approaches. 

Index  Terms — Block  coders,  deblocking,  separable  least  squares, 
variable  projection. 

I.  Introduction 

BLOCK-BASED  coders,  such  as  the  JPEG  standard  for  still 
image  compression  [1]  and  the  MPEG-1  [2]  and  MPEG-2 
[3]  standards  for  coding  of  video  sequences,  are  among  the  most 
common  compression  tools  used  for  compression  of  visual  data. 
Their  low  computational  complexity  along  with  their  good  per¬ 
formance  make  them  a  very  attractive  choice  for  many  applica¬ 
tions.  A  well-known  drawback  of  block  coding  schemes  is  the 
introduction  of  visually  disturbing  blocking  artifacts  in  the  re¬ 
constructed  image.  As  block-based  coders  became  increasingly 
more  popular,  the  need  to  diminish  the  effects  of  blocking  arti¬ 
facts  grew  stronger.  A  great  deal  of  effort  has  been  invested  in 
attempts  to  reduce  these  artifacts,  while  preserving  the  informa¬ 
tion  content  of  the  image. 

Various  postprocessing  techniques  have  been  suggested 
for  the  reduction  of  blocking  artifacts,  based  on  different  ap¬ 
proaches.  Among  these  we  count  techniques  based  on  adaptive 
filtering  [4],  [5],  projection  on  convex  sets  (POCS)  [6]— [8], 
Markov  random  field  modeling  [9],  [10],  and  many  others  (see 
[11]  for  a  review).  Most  of  these  methods  work  well  when  the 
image  is  compressed  at  high  or  medium  bit  rates,  but  the  ma¬ 
jority  of  them  fail  to  produce  satisfactory  results  when  applied 
on  images  compressed  at  low  bit  rates. 

In  search  of  a  more  suitable  scheme  for  low  bit-rate  block 
coders,  several  authors  considered  the  idea  of  downsampling 
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before  coding  and  subsequent  upsampling  at  the  decoder  stage. 
Indeed,  in  [12],  for  instance,  Bruckstein  et  al.  considered  down- 
sampling  an  image  before  applying  the  JPEG  coding  algorithm, 
and  interpolating  at  the  decoder  stage  to  obtain  the  image  in  full 
resolution.  Earlier,  Zeng  and  Venetsanopoulos  [13]  considered 
an  interpolative  image  coding  scheme  based  on  these  notions, 
and  Jung,  Mitra,  and  Mukherjee  in  [14]  introduced  the  idea  of 
’’subband  DCT,”  which  they  later  used  for  image  resizing  in  the 
compressed  domain  [15].  Most  recently,  Dugad  and  Ahuja  [16] 
proposed  an  elegant  scheme  for  spatial  scalability  using  non- 
scalable  coders  that  again  makes  use  of  the  down/up-sampling 
methodology. 

These  approaches  have  several  attractive  properties.  First  and 
foremost,  at  low  bit  rates  there  is  a  marked  gain  in  performance, 
both  in  terms  of  PSNR  and  in  terms  of  visual  quality.  Second, 
the  computational  complexity  involved  in  coding/decoding  is 
substantially  reduced,  since  the  input  to  the  encoder  is  consid¬ 
erably  smaller  in  size.  In  addition,  the  range  of  low  bit  rates 
is  expanded,  allowing  compression  at  lower  bit  rates.  Finally, 
since  most  of  these  methods  do  not  modify  the  basic  coding  al¬ 
gorithm,  but  only  apply  preprocessing/postprocessing,  they  can 
be  used  in  applications  where  the  codec  is  already  implemented 
without  introducing  substantial  modifications. 

Motivated  by  these  features,  the  authors  of  [12]  specifically 
derived  an  analytical  model  for  the  JPEG  encoder  in  order  to 
obtain  an  optimality  criterion  on  the  sampling  factor  for  a  given 
input  image.  Throughout  their  experiments,  they  used  fixed  fil¬ 
ters  for  decimation  and  interpolation,  and  did  not  consider  the 
effect  of  different  filters  on  the  quality  of  the  results.  To  the  best 
of  our  knowledge,  all  of  the  papers  that  have  taken  the  same 
approach  used  fixed,  standard  filters  in  the  down/up-sampling 
process. 

The  main  contribution  of  this  paper  is  to  address  the  issue 
of  finding  optimal  filters  in  practice  for  the  decimation  and  in¬ 
terpolation  stages  in  order  to  achieve  better  performance.  We 
shall  see  that  the  use  of  carefully  chosen  filters,  based  on  least 
squares  (LS)  optimality,  results  in  a  significant  gain  in  perfor¬ 
mance,  both  visually  and  quantitatively.  We  do  not  derive  a  spe¬ 
cific  model  for  the  encoder,  but  consider  it  to  be  a  “black  box.” 
We  view  this  as  an  advantage,  rather  than  a  shortcoming  because 
under  this  framework,  the  results  of  this  paper  are  not  restricted 
solely  to  a  particular  block  coding  scheme,  and  can  be  applied 
to  other  coders  as  well. 

The  rest  of  this  letter  is  organized  as  follows.  Section  II 
presents  the  algorithm  for  finding  the  optimal  interpolation 
filter,  by  formulating  the  problem  as  linear  LS  minimization. 
Section  III  addresses  the  problem  of  jointly  optimizing  both 
filters,  which  we  formulate  as  a  nonlinear,  but  separable  LS 
problem.  We  apply  a  mechanism  known  as  variable  projection 
to  reduce  this  problem  into  a  nonlinear  LS  problem  in  one 
set  of  variables,  and  then  demonstrate  that  even  a  suboptimal 
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Fig.  1.  Sampling/filtering  scheme  for  block  coders. 


approach  results  in  superior  performance.  In  Section  IV  we 
provide  a  discussion  of  the  results,  summarize  our  findings  and 
make  some  concluding  remarks. 


II.  Optimal  Interpolation  by  Linear  Least  Squares 

Throughout  the  rest  of  this  paper,  we  consider  the  system 
shown  in  Fig.  1.  An  input  image  is  convolved  with  a  linear 
filter  f ,  and  then  downsampled  by  a  factor  k.  The  low-resolu¬ 
tion  image  is  then  encoded  using  a  block  coder.  At  the  decoder, 
the  image  is  first  decoded  using  the  block  decoder,  then  upsam¬ 
pled  back  to  its  original  resolution  and  filtered  by  a  filter  g  to 
produce  the  reconstructed  result.  The  authors  in  [12]  took  f  to 
be  a  standard  antialiasing  filter,  and  g  to  be  a  linear  interpola¬ 
tion  kernel  (hat  function).  The  sampling  factor  k  was  chosen 
according  to  their  analytical  predictions.  Throughout  this  paper, 
we  will  assume,  for  simplicity,  k  =  2.  For  a  thorough  discussion 
on  choosing  the  optimal  k,  the  reader  is  referred  to  [12]. 

Let  X  denote  the  input  image  of  size  m  X  n,Y  denote  the 
image  after  decimation  (size  m/2  x  n/ 2),  and  X  the  recon¬ 
structed  image  after  interpolation  as  in  Fig.  1 .  Our  ultimate  goal 
is  to  minimize  the  Z2  error  norm  ||  jf  —  X\\2-  At  first,  we  shall 
only  consider  the  optimization  of  the  interpolation  filter  g,  while 
keeping  f  fixed.  To  do  so,  we  note  that  the  interpolation  stage  can 
be  equivalently  expressed  as  matrix  multiplication.  In  our  for¬ 
mulation  below,  we  consider  the  upsampling  and  filtering  steps 
as  a  unified  process  and  consequently  use  a  set  of  filters,  rather 
than  first  inserting  zeros  and  then  using  one  filter  for  interpola¬ 
tion. 

Specifically,  we  define  X^p,q\p,  q  G  {0, 1},  such  that 


X^\i,j)  =  X(2i+p,2j  +  q) 

1. e.,  XM  is  X  shifted  by  ( p ,  q )  and  then  downsampled  by 

2.  Clearly,  the  set  /X^p'ql  ,p,q  G  {0, 1}}  is  just  a  reordered 
version  of  X.  In  addition,  we  let  x(j,'q>  denote  the  row-stacked 
form  of  X(p,q\  i.e.,  eachx^’9^  has  length  mtt/4,  with  elements 


Now,  for  our  two-dimensional  (2-D)  interpolation  filter  g  with 
dimensions  Z  x  Z,  we  similarly  define  {g(J,’9),0  <p,q  <  l}, 
where  each  g(p’9)  is  a  vector  of  length  Z2  that  represents  the  filter 
in  the  filter  set  which  produces  X^p,q^  by  filtering  Y.  Finally, 


Fig.  2.  Optimal  interpolation  for  Barbara. 


we  construct  a  matrix  $  with  dimensions  mn/A  X  Z2  out  of  the 
image  Y,  of  the  form 


$  = 


<0 

Ki 


4>\ 


(m/2)-l,(n/2)-l 


0) 


where  <pl  •  is  the  row-stacked  form  of  an  Z  X  Z  window,  centered 
around  the  pixel  location  (i,j)  of  Y . 

Using  these  definitions,  we  can  now  express  each  x^p'q'>  (and 
equivalently  X)  as  a  product  of  the  matrix  $  and  the  filter  g(p’q) 
in  vector  form 


(2) 


Looking  at  (2),  we  can  immediately  see  that  an  optimal  solution 
(in  the  LS  sense)  is  obtained  by  setting  X  =  X  and  minimizing 
over  all  g^,g) .  Specifically,  we  solve 


min 

g (P’Q) 


(p,q) 


(3) 


This  is  an  LS  problem,  with  the  solution  g/rqq>  given  by 


g  Ipa)  =  $+x' 


(4) 


where  <1>  denotes  the  pseudo-inverse  of  $  [17],  Obviously, 
in  practice,  there  is  no  need  to  construct  the  matrix  $  nor  its 
pseudo-inverse  .  Rather,  (3)  is  solved  by  applying  recursive 
least  squares  (RLS). 

To  test  the  performance  of  this  method,  we  used  a  standard 
JPEG  coder  as  our  block  coder  and  applied  this  optimization  al¬ 
gorithm  to  several  standard  test  images.  We  also  compared  our 
results  with  those  of  Bruckstein  et  al.  [  12].  Figs.  2-4  show  the 
rate-distortion  curves  obtained  for  the  images  Barbara,  Gold- 
hill,  and  Boats,  respectively.  For  the  decimation  filter,  we  used 
an  1 1-tap  antialiasing  filter,  designed  using  a  Hamming  window 
and  normalized  cutoff  frequency  u>n  =  0.5.  For  the  interpola¬ 
tion  filter,  we  set  1  =  5,  and  solved  (4)  to  find  the  optimal  g(p’q\ 
The  results  clearly  display  a  significant  gain  in  performance  over 
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Fig.  3.  Optimal  interpolation  for  Goldhill. 


Fig.  4.  Optimal  interpolation  for  Boats. 

the  original  results  in  [12].  We  see  that  the  optimal  curve  inter¬ 
sects  the  JPEG  curve  at  a  higher  bit  rate  than  the  curve  obtained 
by  Bruckstein  et  al.  This  essentially  means  that  our  algorithm  is 
applicable  to  a  wider  range  of  bit  rates,  since  it  performs  better 
than  direct  JPEG  compression  up  to  a  higher  bit  rate.  The  visual 
improvement  over  [12]  is  clearly  evident  in  Figs.  5-7.  We  can 
see  that  the  optimal  filter  provides  a  significantly  sharper  image, 
while  virtually  eliminating  the  blockiness.  In  terms  of  PSNR,  in 
all  three  cases,  when  coded  with  0.2  bpp,  the  gain  over  standard 
JPEG  is  around  1-1.5  dB. 

It  would  be  interesting  to  examine  the  filters  obtained  in  the 
optimization  process.  Fig.  8  shows  the  2-D  frequency  response 
of  the  optimal  interpolation  filters  obtained  for  Barbara  (recall 
that  4  filters  are  needed).  As  we  can  see,  all  four  filters  exhibit 
a  low-pass  behavior  with  some  high-pass  components.  In  other 
words,  the  filters  attempt  to  smooth  out  the  blockiness  while 
preserving  edge  and  texture  information. 

The  overhead  of  sending  the  filter  coefficients  to  the  decoder 
is  not  included  in  the  rate  calculation,  yet  it  is  in  the  order  of 


(c) 


Fig.  5.  Compression  results  for  Barbara ,  0.2  bpp.  (a)  JPEG,  PSNR  = 
23.42  dB.  (b)  Bruckstein  et  al.,  PSNR  =  24.19  dB.  (c)  Optimal  interpolation, 
PSNR  =  24.74  dB. 

a  few  hundred  bits  per  filter,  which  is  negligible  even  for  very 
low  bit  rates  (we  found  that  for  the  parameters  used  in  our  ex¬ 
periments,  each  filter  would  occupy  around  200  bits  when  prop¬ 
erly  entropy  coded).  We  would  like  to  comment  that  the  4-fil- 
ters  approach  we  have  taken  is  completely  equivalent  to  a  single 
filter  approach,  in  which  case  zero-padding  is  used  prior  to  fil¬ 
tering,  and  only  one  vector  x  is  constructed.  Indeed,  it  is  not 
hard  to  see  that  the  LS  problem  4>g  =  x,  with  $  a  row-stacked 


IEEE  TRANSACTIONS  ON  CIRCUITS  AND  SYSTEMS  FOR  VIDEO  TECHNOLOGY,  VOL.  15,  NO.  1,  JANUARY  2005 


157 


Fig.  6.  Compression  results  for  Goldhill,  0.2  bpp.  (a)  JPEG,  PS  NR  = 
27.43  dB.  (b)  Bruckstein  etal.,  PSNR  =  27.95  dB.  (c)  Optimal  interpolation, 
PSNR  =  28.91  dB. 

form  of  the  image  Y  after  zero-padding,  can  be  decoupled  into 
four  subproblems  of  the  form  of  (3).  The  advantage  of  using  the 
4-filters  approach  is  apparent  when  the  matrix  $  needs  be  con¬ 
structed  explicitly,  in  which  case  it  is  four  times  smaller  than  if 
zero-padding  was  used.  However,  in  instances  where  the  sam¬ 
pling  factor  is  greater  than  two,  or  when  it  is  a  rational  number, 
we  note  that  such  an  approach  may  not  be  feasible,  and  so  the 
problem  should  be  formulated  using  a  single  filter. 


Fig.  7.  Compression  results  for  Boats ,  0.2bpp.  (a)  JPEG,  PSNR  =  2G.S2  dB. 
(b)  Bruckstein  et  at.,  PSNR  =  27.42  dB.  (c)  Optimal  interpolation,  PSNR  = 
28.58  dB. 

III.  Optimal  Decimation  and  Interpolation  by  Variable 
Projection 

Inspired  by  the  results  obtained  when  optimizing  over  the  in¬ 
terpolation  filter,  we  now  turn  our  attention  to  the  decimation 
filter.  If  we  consider  minimizing  the  difference  between  X  and 
X  over  both  f  and  g.  then  the  image  Y  at  the  output  of  the  de¬ 
coder  is  no  longer  fixed,  hence  our  matrix  #  of  (1)  which  origi- 
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Fig.  8.  Optimal  interpolation  filters  obtained  for  Barbara. 

nated  from  Y  is  now  dependent  on  f,  i.e.,  $  =  4>(f).  Nonethe¬ 
less,  we  can  still  write  our  problem  as  a  matrix  system  in  the 
form  of  (2) 

X^’9)  =  *(f)gfr-9>  (5) 

and  our  minimization  problem  (3)  now  becomes 


min  y^P^)  _  $(f)g(p’9)||2-  (6) 

ftg(p,«)  2 

This  is  a  nonlinear  LS  problem  with  respect  to  the  variables 
f ,  g(P’9) .  To  find  a  solution  for  this  problem,  we  apply  Golub 
and  Pereyra’s  Variable  Projection  method  [18J.  This  method  is 
applicable  when  dealing  with  separable  nonlinear  LS  problems, 
i.e.,  problems  for  which  the  model  function  is  a  linear  combi¬ 
nation  of  nonlinear  functions.  Since  first  published,  the  variable 
projection  method  (sometimes  coined  VARPRO)  has  been  ap¬ 
plied  in  a  variety  of  different  fields  in  engineering  and  scientific 
computing  (see  [19]  for  an  extensive  survey). 

In  our  formulation,  we  have  two  separable  sets  of  variables, 
namely  f  and  g Moreover,  the  dependence  on  g(,v/)  is 
linear.  Hence,  we  can  use  VARPRO  to  transform  the  problem 
into  a  nonlinear  LS  problem  in  one  set  of  variables.  More  specif¬ 
ically,  assume  for  the  moment  we  know  the  optimal  f .  If  we  plug 
this  f  into  (5),  then  <&(f)  is  now  fixed,  hence  we  again  face  a 
linear  LS  problem,  with  the  solution  readily  given  by 


=$(f)+x(^).  (7) 


Frequency  -1  -1  Frequency 


Now,  we  can  use  this  expression  for  gfr'i)  in  our  minimization 
problem  (6),  leading  to 

min  -  $(f)$(f)+x^  ~  =  mm  P^(f)x^’«)  * 

(8) 

where  P^(f)  =  I  —  $(f)$(f)+  is  the  projector  on  the  orthog¬ 
onal  complement  of  the  column  space  of  3>(f). 

As  we  can  see,  by  using  VARPRO,  we  have  essentially 
eliminated  the  minimization  with  respect  to  g^p,q\  and  we  are 
left  with  a  nonlinear  LS  problem  with  respect  to  the  decimation 
filter  f .  Evidently,  this  procedure  not  only  reduces  the  dimension 
of  the  parameter  space  but  also  results  in  a  better-conditioned 
problem  (see  [19]  for  details).  Nonetheless,  solving  (8)  is 
still  a  difficult  task,  due  to  the  nonlinearity  of  the  problem 
at  hand.  For  our  purposes,  we  shall  pursue  a  near-optimal 
approach,  by  restricting  our  discussion  to  a  fraction  of  the 
parameter  space  for  this  optimization  problem.  We  notice  that 
the  parameter  space  for  this  problem  is  /-dimensional,  where  / 
is  the  length  of  the  filter  f  (assuming  separability  of  f).  Rather 
than  considering  the  entire  space,  we  consider  a  family  of 
lowpass  filters  with  a  varying  cutoff  frequency.  For  the  design 
of  the  lowpass  filter  fr p(cj).  we  use  the  window  method  with 
a  Hamming  window  to  design  a  separable  filter  with  cutoff 
frequency  u.  Consequently,  a  suboptimal  solution  of  (5)  is 
given  by  the  lowpass  filter  fr p(u>),  where  Cj  is  given  by 

u>  =  argmm  P i(fLpM)x(p’9)  ^  (9) 

This  is  a  univariate  minimization  problem,  which  we  solve  using 
Brent’s  classic  algorithm  for  minimization  without  derivatives 
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Fig.  9.  Optimal  interpolation  and  decimation  for  Barbara. 


Fig.  10.  Optimal  decimation  and  interpolation  for  Barbara ,  0.2  bpp,  PS  NR  = 
25.5  dB. 

[20].  We  note  that  the  optimization  problem  (9)  is  solved  once 
at  the  encoding  stage,  and  the  resulting  decimation  filter  need 
not  be  sent  to  the  decoder  as  it  is  not  used  in  the  decoding 
process.  Hence,  no  additional  overhead  is  introduced  by  this 
stage. 

Before  examining  the  experimental  results,  we  shall  comment 
about  the  computational  burden  involved  in  each  iteration  of 
Brent’ s  algorithm.  Each  such  iteration  requires  coding  the  image 
after  applying  the  filter  fL p(u>),  followed  by  a  decoding  stage 
and  computation  of  the  residual  in  (9).  In  our  experiments  on 
typical  512  x  512  images,  the  execution  time  per  iteration  was 
3.92  s,  when  executed  on  a  P4  2-GHz  machine,  using  Matlab.  As 
mentioned  above,  the  solution  of  (9)  is  carried  out  once  during 
the  encoding  stage,  and  no  further  computations  need  be  done 
at  the  decoder. 

Fig.  9  shows  the  rate-distortion  curve  obtained  when  applying 
this  algorithm  on  the  test  image  Barbara.  Clearly,  optimizing 
over  both  filters  results  in  significantly  better  quality  compared 


to  optimization  of  the  interpolation  filter  alone.  This  is  also  in¬ 
herent  in  Fig.  10,  which  displays  the  visual  results  for  a  bit  rate 
of  0.2  bpp  (compare  Fig.  5).  The  image  obtained  from  the  joint 
optimization  is  sharper  and  exhibits  more  details,  and  is  free  of 
any  blocking  artifacts.  Interestingly  enough,  the  optimal  deci¬ 
mation  filter  found  in  this  case  is  a  low-pass  filter  with  cutoff 
frequency  ui  =  0.97,  which  is  essentially  an  identity  filter.  In¬ 
tuitively,  this  means  that  we  do  not  want  any  filtering  prior  to 
downsampling,  in  order  to  preserve  the  texture  that  dominates 
Barbara. 

IV.  Discussion  and  Conclusion 

This  letter  presented  an  algorithm  for  finding  optimal  filters 
in  a  sampling/compression  scheme.  By  a  LS  argument,  we  have 
shown  that  the  use  of  optimal  interpolation  filters  achieves  a 
significant  gain  in  performance,  compared  to  using  common  fil¬ 
ters.  We  demonstrated  that  a  near  optimal  approach  for  finding 
the  decimation  filter  also  results  in  a  marked  improvement  in 
quality.  In  addition,  by  using  optimal  filters,  we  outperform 
JPEG  for  a  wider  range  of  bit  rates  compared  to  [12],  making 
this  sampling/compression  scheme  applicable  in  more  diverse 
situations. 

We  note  that  even  greater  improvements  can  be  gained  by  tai¬ 
loring  the  filter  optimization  to  the  block  structure  of  the  coder. 
By  using  one  set  of  filters  for  block  boundaries,  and  another  for 
block  interior,  we  are  likely  to  get  better  performance,  as  there 
is  much  to  be  gained  by  treating  these  two  regions  separately. 
Furthermore,  we  note  that  the  restriction  of  the  parameter  space 
when  solving  (9)  to  separable  low-pass  filters  is  quite  stringent. 
Better  performance  may  result  by  considering  a  wider  variety 
of  decimation  filters,  while  keeping  the  optimization  problem 
tractable.  These  issues  are  left  for  future  research. 
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